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Abstract 

We review and extend results for mutation, selection, genetic drift, and migration in a one- 
dimensional continuous population. The population is described by a continuous limit of the step- 
ping stone model, which leads to the stochastic Fisher-Kolmogorov-Petrovsky-Piscounov equation 
with additional terms describing mutations. Although the stepping stone model was first proposed 
for population genetics, it is closely related to "voter models" of interest in nonequilibrium sta- 
tistical mechanics. The stepping stone model can also be regarded as an approximation to the 
dynamics of a thin layer of actively growing pioneers at the frontier of a colony of microorganisms 
undergoing a range expansion on a Petri dish. We find that the population tends to segregate into 
monoallelic domains. This segregation slows down genetic drift and selection because these two 
evolutionary forces can only act at the boundaries between the domains; the effects of mutation, 
however, are not significantly affected by the segregation. Although fixation in the neutral well- 
mixed (or "zero dimensional") model occurs exponentially in time, it occurs only algebraically fast 
in the one-dimensional model. If selection is weak, selective sweeps occur exponentially fast in 
both well-mixed and one-dimensional populations, but the time constants are different. We also 
find an unusual sublinear increase in the variance of the spatially averaged allele frequency with 
time. Although we focus on two alleles or variants, g-allele Potts-like models of gene segregation 
are considered as well. Our analytical results are checked with simulations, and could be tested 
against recent spatial experiments on range expansions off linear inoculations of Escherichia coli 
and Saccharomyces cerevisiae. 

PACS numbers: 87.23.Kg, 87.23. Cc, 87.18.Hf, 64.60.De 

Keywords: stepping stone model, stochastic Fisher-Kolmogorov-Petrovsky-Piscounov equation, selective 
sweep, voter model, Eden model. 
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I. INTRODUCTION 



The quantitative theory of evolution is an important open problem. The theory of evo- 
lutionary dynamics is necessary to determine the history of species migrations, and it could 
shed light on the origin and development of hfe. Moreover, a better understanding of the 
evolutionary dynamics could help control epidemics, fight diseases with an evolutionary char- 
acter such as cancer and acquired immune deficiency syndrome, and guide the engineering 
of artificial evolution for practical applications. 

Most of the current understanding of evolutionary dynamics comes from population ge- 
netics, a scientific discipline that studies how evolutionary forces shape the genetic diversity 
of populations. The majority of theoretical models and experiments in population genetics 
study only one or a few well-mixed populations, that is populations without spatial struc- 
ture, where every individual is equally likely to interact with any other individual inside 
the same population. Microorganisms growing and evolving in a well-mixed liquid culture 
provide an important example. While nonspatial models are often easier to analyze than 
spatial ones, they do miss what can be essential features of natural populations. 

In nature, organisms often occupy areas that are much larger than the square of the disper- 
sal distance, that is the distance typically traveled by an individual in one generation. This 
causes two main problems for well-mixed-population models. First, well-mixed-population 
models underestimate the role of genetic drift (fluctuations due to the discreteness of the 
number of individuals). The difference arises because the organisms can only interact with 
their neighbors, and the number of neighbors within the dispersal distance is much smaller 
than the total number of organisms in the entire population. Second, well-mixed-population 
models neglect the spatial structure of the population that can be created by external factors 
or by internal dynamics. Such spatial structures often exist, and, as we show in this paper, 
they can significantly affect evolutionary processes in the population. 

Well-mixed-population models are particularly inadequate when applied to expanding 
populations. Expansions are very common in biology. Species spread to new territories 
from the locations where they first evolved. Expansions also occur because of environmental 
changes such as the global warming and the glacial cycles or due to sudden long distance 
migrations to new habitats. Even though well-mixed-population models can account for the 
growing number of individuals (population size), these models do not capture the fact that 
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the newly settled areas are colonized by the offspring of only a small number of individuals 
at the expanding front. Since the ancestral population is small, the genetic drift is strong. 
As a result, neutral genetic diversity decreases with the distance from the origin of the 
expansion. This reduction in genetic diversity, which is often called "the founder effect" l| , 
has been observed in humans [2, Q] and many other species. For example, the founder effect 
in the population waves following the receding glaciers is believed to be responsible for the 
reduced genetic diversity in high latitude regions compared to equatorial ones |4|. 

The spreading of Escherichia coli (E. coli) and Saccharomyces cerevisiae (S. cerevisiae) 
on Petri dishes has been investigated in recent experiments by Hallatschek et al. [5]. In 
these experiments, microbes grown in the dark carried one of two selectively neutral alle- 
les, differing only in a gene encoding for proteins with two distinct fluorescence spectra. 
Figure [1] shows the expansion of an initially well-mixed 50 — 50 population of E. coli into 
two unoccupied half planes initiated by a razor blade inoculation with cells grown up in 
liquid culture. The distinctive feature illustrated by the typical experiment in Fig. [T] is that 
the population does not remain well-mixed; instead, it segregates into well-defined domains. 
The segregation occurs because the strong genetic drift associated with reduced population 
size facilitates fixation of one of the two alleles at the front. 

Analogous phenomena should also occur in a nonexpanding one-dimensional population 
because its dynamics is similar to the dynamics of the front of a growing population. The 
front of a population wave and a literally one- dimensional habitat are not exactly equivalent 
because the contour of the front undergoes undulations while a one- dimensional habitat has 
a fixed linear shape. Nevertheless, both are effectively one-dimensional and should deviate 
from the predictions of well-mixed-population models in similar ways. The advantage of 
a fiat one-dimensional habitat is that it is easier to analyze. In addition, although most 
species live in effectively two dimensional habitats, a quasi one-dimensional habitat could 
describe a bank of a river, a sea coast, and a slope of a linear mountain range. 

To study the dynamics of a population analytically, we adopt the stepping stone model 
proposed by Kimura and Weiss Gj. This model considers many well- mixed populations, 
demes, located on a spatial lattice. Each deme is subject to mutation, selection, genetic 
drift, and short range migration between neighboring demes. In the limit of weak evo- 
lutionary forces and large number of demes, the stepping stone model is equivalent to the 



continuous models proposed by Wright \^ and Malecot 



and is described by the stochastic 
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FIG. 1: (Color online) Spatial segregation in an expanding microbial population. Different colors 
label different alleles. Initially well-mixed population occupies a narrow horizontal linear region 
lying between the arrows, which show the direction of the growth. As this population expands, 
it segregates into well defined monoallelic domains. Details of the experiment are presented in 
Ref. [5|. 

Fisher-Kolmogorov-Petrovsky-Piscounov equation [o], [l^ with additional terms representing 
mutation. On the other hand, when each deme contains only one organism, the model is 
analogous to the Eden model ll| used to describe the growth of interfaces and the voter 
model Q. 

We also performed numerical simulations to better understand the relationship between 
the experiments in Ref. 5| and our analytical results. An illustrative simulation (with 
periodic boundary conditions) is shown in Fig. O which also shows the difference between a 
growing population front with undulations and a literally one- dimensional habitat advancing 
uniformly in time. Figure [3] shows qualitative agreement between the experiments and the 
simplified row-by-row growth model that we studied analytically. 

In this paper, we first focus on the spatial segregation due to genetic drift and its effect 
on the dynamics of a one-dimensional population. We find that segregation of two neutral 
alleles has two stages. During the first stage, distinguishable domains emerge from the well- 
mixed population. During the second stage, domain boundaries diffuse and annihilate upon 
collision. As a result, some of the domains vanish whereas others grow. We show how our 
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FIG. 2: (Color online) Computer simulations. The blank hexagons represent empty sites, and 

different colors of the occupied hexagons represent different alleles, (a) The model of an expanding 

population front. The highlighted hexagon is a randomly chosen cell that can reproduce and deposit 

an identical offspring in any of its four empty nearest neighbor sites (shown with arrows) with equal 

probability, (b) The model of a one-dimensional habitat. Each row represents a generation so each 

row is completed before moving on to the next one. An empty site can be filled only by an offspring 

of one of its nearest neighbors in the previous generation (shown with arrows). Both (a) and (b) 

show the effects of genetic drift (sampling error) when, e.g., the second from the left cell in the 

bottom row leaves no offspring. Such events lead to coarsening seen in (c) and (d). (c) and (d) 

are single simulation runs for models in (a) and (b) respectively. A population of 100 cells was 

wrapped around a cylinder to illustrate periodic boundary conditions used in this paper. Note that 

in (d) the front is flat whereas in (c) it is rough. This roughness affects some aspects of the shapes 

of the monoallelic domains shown in (c): A domain boundary followed from its lowest point to its 

highest point always goes up in (d), but, in (c), it sometimes goes down. As discussed in Ref. Q], 

domain walls are expected to wander more vigorously in (c) than in (d). 
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FIG. 3: (Color online) Qualitative comparison of a gene segregation experiment from a linear 
inoculation (inset) and the simulation of a one-dimensional habitat. 



calculations might be used to extract the diffusion constant and the effective population size 
from experiments like those in Ref. and discuss how well the model describes the behavior 
of microbes. A detailed comparison (beyond the qualitative agreement we find with the main 
features) would require more extensive and precise experiments; we hope such experiments 
will be carried out in the future. The spatial segregation dramatically changes the effects 
of genetic drift and selection on the population compared to the predictions of well-mixed- 
population models. For the neutral model without mutation, we find that local diversity 
or "heterozygosity" decays as and the standard deviation of the global fraction of 

an allele grows sub diffusively as t^^^. We also study the dynamics in the presence of weak 
selection and find that it differs markedly from that of a well-mixed population. Because 
of the spatial segregation into domains, selection acts only near domain boundaries, which 



constitute only a small fraction of the population. Hence, extinction of a deleterious allele 
proceeds much more slowly in one-dimensional populations than in well-mixed populations. 
Unlike genetic drift and selection, the effects of mutation in the spatial model are essentially 
the same as in the well-mixed-population model, but the spatial model gives a more accurate 
description of the population and accounts for the spatial correlations. 

A substantial fraction of our results for the neutral dynamics in a one-dimensional habitat 
has been derived previously in population genetics 6, Q, Q, ecology 



librium statistical mechanics 



12, 



15], and nonequi- 



16(1 . Here, we present a single self-contained derivation of 



these earlier results in a novel context of expanding populations in two dimensions and in a 
language familiar to physicists, with future microbial tests of the theory in mind. Our new 
results are primarily confined to the analysis of natural selection. 

This paper is organized as follows. First, we review classical results for well-mixed pop- 
ulations in Sec. [Tll We then introduce the one- dimensional stepping stone model in Sec. IIIII 
and derive the equations of motion for spatial correlation functions. In Sec. HV] and Sec. |V] 
we solve these equations for zero and nonzero mutation rates respectively. While the neutral 
stepping stone model has been treated before, we derive some new results and use a differ- 
ent technique that can be easily extended to radially expanding populations. The effects 
of selection are considered in Sec. IVTl and in Sec. IVIII we test our analytical results with 
simulations. Various details are relegated to Appendices RllDl In Appendix |El we indicate 
how some of the 2-state (i.e., "2-allele") results can be generalized for the Potts-model-like 
nonequilibrium dynamics of g-alleles with g > 3. 



II. POPULATION GENETICS IN WELL-MIXED POPULATIONS 

Well-mixed-population models are relevant to microorganisms vigorously shaken in a test 
tube, but they do not describe spatial phenomena. Indeed, if cells visit all parts of the test 
tube during a cell division time, they are effectively zero-dimensional. Nevertheless, they 
can serve as a useful reference point to which spatial models can be compared. Well-mixed 
populations also provide a simple context to introduce genetic drift, mutation, and selection; 
and the stepping stone model presented in Sec. IIIII uses a well-mixed-population model to 
describe the dynamics of allele frequencies within the demes. This section summarizes the 
classical results of nonspatial population genetics, which are primarily due to Wright, Fisher, 
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Haldane, and Kimura; Refs. [17| and [l8j] provide a good introduction to the subject and 
refer to the original hterature, which is too extensive to be discussed here. 

To simphfy the discussion and to make a direct connection with the experiments in 
Ref. j^, we consider two alleles in a population of haploid organisms. (For N diploid 
organisms, the theory is essentially the same under certain assumptions, provided one fo- 



cuses on the dynamics of 2N gene copies in each generation; see Ref. [17|.) The two-allele 
approximation may seem very restrictive, but many of our results can be generalized to an 
arbitrary integer number g > 3 of alleles. In addition, a two-allele model can be used to 
describe the dynamics of an allele of interest (with or without a selective advantage) when 
all other alleles are neutral. We assume that each of the individuals in the population can 
die, give birth (divide), and mutate. The details of this birth and death process are species 
dependent, but the dynamics on time scales larger than the generation time Tg is believed 
to be universal provided N is large. This universal dynamics is often referred to as the 
diffusion or continuous approximation. Two simple models are commonly used to illustrate 
the continuous approximation: the Wright-Fisher model and the Moran model. Here, we 
use the latter because it more closely resembles microbes with overlapping generations. 

First, we consider the Moran model without selection and mutation. During a time 
step, two individuals are randomly selected with replacement from the population. The 
first individual is chosen to reproduce, and the second one to die; thus the total number 
of the organisms is conserved. If the "frequency" of allele one (i.e., the fractional number 
of individuals with genotype one) at time step t is fit), then, at the next time step, it 
is / + with probability /(I — /), / — with probability /(I — /), and / with 
probability + (1 — /)^. The expectation value and variance of f{t + 1) are then given by, 

(/(t + l)) = /(t), (1) 

{[f{i+ 1) - (/(t + 1))]^) = (2) 

where angular brackets represent average with respect to the random choice of individuals 
for reproduction and death. Because only one of organisms gives birth in a Moran time 
step, t measures time in fractional generation time, Tg/N. 

Equations ([1]) and ([2]) imply that f{t) performs an unbiased random walk in the space 
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of allele frequencies. In the continuum limit, this random walk can be described by the 

nn 

following Fokker-Planck equation with a frequency dependent diffusion coefficient |17l. |18| 

^^ = ^|j[/(l-/W,/)], (3) 

where P{t, /) is the probability density function for / at time t measured in generations, 
and Dg is the genetic diffusion constant. Here, t is the time measured in generations; as 
discussed above, N Moran time steps constitute a generation time Tg. Thus, in the Moran 
model, we have 



lead to an iden- 



Alternative reproduction schemes, such as Wright-Fisher sampling, [17|, 
tical equation, with a different numerical coefficient in Eq. (jl]). 

Equation ([3]) is subject to absorbing boundary conditions ^ at / = and / = 1 because, 
if one of the alleles is lost, it cannot appear again in the absence of mutation. Therefore, 
the population eventually becomes fixed at one of the absorbing states. We calculate the 
rate of the fixation by considering the average heterozygosity of the population 

H{t)^{h{t)) = {2m[i-fm, (5) 

which is the (averaged over realizations) probability that two randomly selected individuals 
have different alleles. When the population is close to the fixation (/ ^ or / ^ 1), the 
heterozygosity is close to zero. The equations of motion for F{t) = (fit)) and H{t) follow 
from Eq. ([3]) by multiplying both sides with / or h, integrating over /, and eliminating the 
derivatives with respect to / via integration by parts. The results are 

=0, (6) 



dt 

dH{t) _ 
dt ~ 



-^gHit). (7) 



^ Since Eq.[3]is singular at the boundaries, we require lim/_^ooO,i /(I — f )P{^-, f) — 0- See Refs. p^] and Q 
for a more detailed discussion. 
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Equations ({6]) and ([7]) imply that, while the average frequencies of these neutral alleles 
do not change F = {f) = f(t = 0) = Fq, the population reaches fixation exponentially 
fast, H{t) = H{0)e-^^' = Fo(l - Fo)e-^«*. 

The average heterozygosity is closely related to the variance of f{t), the fraction of the 
first allele, 

v{t) = nm - {f{tm = Fim - m] - Imt). (s) 

Thus, even if a population starts with zero variance, the fluctuations grow until the variance 
reaches its maximum value of -Fo(l ~ -fo), which corresponds to a population fixed to allele 
one with probability Fq and to allele two with probability 1 — Fq. Note that, for small t, 
V{t) grows linearly with time, but, at large times, the variance approaches its limiting 
value exponentially fast. The linear growth of variance at small times also follows from the 
Fokker-Planck equation because, at small times, Eq. ([3]) can be approximated by a diffusion 
equation with constant diffusivity. 

Next, we generalize Eq. ([3]) to account for mutations. In the Moran model, mutation is 
included at the end of a time step by allowing the offspring to mutate with probability /ii2 
from allele one to allele two and with probability /i2i from allele two to allele one. If the 
frequency of allele one at time step i is f{i), then, at the next time step, the expectation 
value of /(t + 1) is given by 

(/(f+l)> = /(f) + ^21^^^<|^^. (9) 

and the variance of /(t + 1) is given by Eq. ([2]) to the leading order in the mutation rates 
and the inverse population size. 

Since the expectation value of f{t) changes with time, mutation leads to an /-dependent 
drift term in the Fokker-Planck equation. Upon recalling that N Moran time steps equal 
one generation time, we have 



^' - ^ {K - + fl!l)/|P((, /)} 



dt df 



+ ^^[fa-f)PitJ)] 
where /ii2 = p'uTg^ and /ii2 = p'2i'Tg^ are the mutation rates per generation 
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Because the alleles can mutate into each other, the probability flux through the boundaries 
must be zero, so Eq. ffTOj) has reflecting boundary conditions, and a nontrivial stationary 



solution for P{t, /) exists. While the stationary distribution can be obtained easily [18|], it 
is sufficient to analyze the moments F{t) and H{t) introduced above. This will also allow 
us to make a direct comparison with the corresponding solutions of the one-dimensional 
stepping stone model in Sec. |Vl The equations of motion for F(t) and H(t) are obtained 
from Eq. ffTOj) in the same way as for the absence of mutation. The results are 

^^=/i21-(/il2 + /i2l)F(t), (11) 
= -i^g + 2/ii2 + 2l22l)H{t) + 2[^21 + (/il2 - l^2l)F{t)]. (12) 

Since these equations are linear differential equations with constant coefficients, the equilib- 
rium is approached exponentially fast. The stationary solutions, which are obtained in the 
limit t oo, are given below 

F(oo) = -J^, (13) 

/^12 + /^21 

^(oo) ^ ^^<°°>" (14) 

~'~ 2(^ti2+A12l) 

From Eqs. f|T^ and ([8]), we see that, when the population size is large enough, i.e. Dg <^ 
(/ii2 + At2i)) H{oo) ^ 2F(oo)[l — F(oo)], the stationary value of the heterozygosity is con- 
sistent with f{t) ^ F(oo). Thus V{oo) ~ 0, and the fluctuations of f{t) are negligible. 
In the opposite limit, H{oo) oc ^^^^^^^ is signiflcantly smaller, which suggests that most of 
the time the population is flxed to one of the alleles, and mutations lead to rare transitions 
between states with / = and / = 1. Our interpretation of Eq. (HM is consistent with a 



more rigorous analytical and numerical analysis by Duty |21 |. 

Finally, we introduce Darwinian natural selection, which is usually related to the differ- 
ence in the reproduction or survival probability of the organisms. In the continuous time 
limit considered here, both mechanisms of selection lead to the same dynamics; therefore, 
we only consider selection due to different growth rates. In the Moran model, a growth 
rate difference is embodied in modifled probabilities of reproduction: the individual with 
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allele one is chosen to reproduce not with probability / but with probability , 
where wi and W2 are the fitnesses (i.e., growth rates) of alleles one and two respectively. In 
the absence of mutations, this modification results in 

N{wif{t) + W2[l - f{t)]} 
When selection is weak, that is \wi — W2\ <^ wi + W2, Eq. ( |T5l) reduces to 

(/(t + l)) = /(t) + ^/(t)[l-/(t)], (16) 

where s = 2{wi — W2) / {1^1 + 102) is the selective advantage of allele one, which has to be much 
smaller than one for the approximation to hold. When s > 0, allele one is advantageous; 
for 5 < 0, it is deleterious. In the following, we assume that allele one is advantageous 
because one can always relabel the alleles to satisfy this condition. 

Similar to the case of mutations without selection, the variance of f{t + 1) is given by 
Eq. ([2]) to the leading order in s and A^~^, and the Fokker-Planck equation at the level of 
the generation time Tg acquires an /-dependent drift term due to selection: 



+ ^gJ-2[f{^-f)P{t,f)], 

where s = st~^. The equations for F and H are not as useful as before because the hierarchy 
of the moment equations does not close. Nevertheless, Eq. 0171) can be easily analyzed in 
two limits. When the population size is large {Dg ^ s), fiuctuations are not important, 
and ^ ^ sF{t)[l - F{t)]. Upon setting Fq = F(0), we have 



so the selective sweep is exponentially fast. When the fiuctuations dominate the dynamics, 
the selection slightly increases the odds of fixation of the advantageous allele, but does not 
significantly affect the rate of fixation. For a detailed analysis of Eq. ( IT71) see Ref . 13] • 

In the continuous limit, the population genetics of a well- mixed effectively zero- 
dimensional population with genetic drift, selection, and mutation is summarized by the 
following Fokker-Planck (or forward Kolmogorov) equation: 
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-^{[/^21-(/il2+/i2l)/]P(t,/)} (19) 



Although this formulation is appropriate for nonspatial models, an alternative formulation 
via a stochastic differential equation can be generalized to the spatial models more easily. 
Equation (fT9|) is equivalent to 



dm 

dt 



Sf{t)[l - fit)] + /i21 - (^12 + 1^21) fit) + \^gfit)[l - fit)]Tit) ilto), (20) 



(r(ti)r(t2)) = 5(tl-t2), (21) 

where r(t) is a white, zero mean Gaussian noise, and 5(t) is Dirac's delta-function; to get 
the correct Fokker-Planck Equation (fT9l) . one must use Ito's prescription to define how 
Eq. ( l20i) steps the dynamics forward in time. This interpretation of the noise term ensures 
that fit) depends only on r(t') with t' < t as it is appropriate for population genetics. Ito's 
prescription is adopted throughout t he p aper, and a brief introduction to the Ito calculus is 
given in Appendix [D] (see also Refs. [l9|, l2l|, [2^). In Sec. IIII[ we use Eq. f l20|) to formulate 



the stepping stone model in one dimension. 

Well-mixed-population-models do not describe migration and subdivision of natural pop- 
ulations \x\- To remedy this deficiency, two common approaches exist: to assume a uni- 
formly populated spatial habitat with free diffusion or to assume a patchy habitat with a 
prescribed pattern of limited migration between the patches. The former is the subject of 
this paper, and can be regarded as the continuum limit of the stepping stone model 6|, see 



Sec. mil The simplest variant of the latter approach is known as the island model 23|] . The 
island model assumes that all patches or islands have the same number of organisms and 
populations in every patch obey well-mixed-population dynamics. The migration occurs 
between any two patches with equal probability, so, in some sense, this is a mean field or 
infinite-dimensional model. The island model successfully predicts that the organisms are 
more likely to be related locally than globally, but most of its predictions are similar to 
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those of well-mixed-population models because the migration does not account for spatial 
structure. In the limit of an infinitely large number of islands, the effect of migration in and 
out of any patch is equivalent to an effective mutation rate; however, this is not the case in 
a one-dimensional model considered below. 

III. ONE-DIMENSIONAL STEPPING STONE MODEL 

In Sec. mi we formulated a model to describe genetic drift, mutation, and selection in 
an effectively zero- dimensional habitat. For a one- dimensional population considered in this 
section, we extend the model to account for short range migrations during every generation. 
Migration is usually modeled either as exchange of individuals between neighboring island 
populatp,. Memes) 0. Q or as d.persal of offsprmg or adults w,th„r a coufrruous popula- 
tion [7|, Il3|, |2J] . Although the first approach was developed to model patchy populations, it 
can be used to describe continuous populations if the deme sizes are much smaller than the 
whole population, and spatial variations are gradual. In this limit, both migration models 
should give essentially the same results. Here, we start with first approach because it is 
conceptually simpler. 

To specify the one-dimensional stepping stone model, we consider an infinite set of demes 
arranged on a line. Neighboring demes are separated by distance a and indexed by an 
integer I = — oo,...,— I,0,l,...,cx3. Each deme has organisms (but the total population 
size is infinite), and the frequency of allele one in deme / is fi{t). Migration occurs only 
between nearest neighbors, and, every generation, a deme exchanges mN individuals with 
its right neighbor and rfiN individuals with its left neighbor. We assume that the exchange 
fraction m is much smaller than one, and that the individuals of both allelic types are equally 
likely to be exchanged. Thus, in one generation, (/;) changes by rh{fi_i + — 2fi) due 
to migration. The variance of /; grows due to randomness in the exchange process, but this 
increase is negligible compared to the genetic drift within an island. In the continuous time 
limit, fi(t) obeys the following generalization of Eq. (120|) : 



^ = m(/^„i + - 2fi) + sfiil - fi) +/i2i - (/ii2 + /i2i)/i + \/^gfi{l - fi)ri (ltd), (22) 
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{Ti,ih)Ti,it2)) = 5i,i,5{t,~t2), (23) 

where m = mT~^ and 6ij^i^ is Kronecker's delta. We can also write Eq. ( !22l) in the continuous 
space limit by introducing a spatial coordinate x = la, 

2 . 



df ^ d'f 



+ sf{l - /) + - + /i2i)/ + \ DJ{1 - f)T (ltd), (24) 



dt dx 



(r(ti, xi)T{t2, X2)) = 5iti - t2)5{xi - X2), (25) 

where the spatial and genetic diffusion constants are Dg = ma? and Dg = aDg = 
respectively. Thus, the continuous time and space limit the stepping stone model is described 
by the stochastic Fisher-Kolmogorov-Petrovsky-Piscounov equation j^, [l^ with additional 
terms describing mutation. 

Similar to the analysis of the well-mixed-population model discussed in Sec. [Tll we use 
equal-time correlation functions of f{t, x) to characterize the dynamics of the stepping stone 
model. The spatial versions of the average frequency and heterozygosity are defined as 
follows: 

F(t,x) = (/(t,x)), (26) 

H{t,Xi,X2) = (/(t,xi)[l - /(t,X2)]) + {f{t,X2)[l - f{t,xi)]). (27) 

The equation of motion for F(t, x) depends on H{t, x, x), and is readily derived by averaging 
Eq. (HM . which gives 

OF d^F s 

— = Ds-^ + /i2i - (Aii2 + ^l2l)F + -H{t, X, x). (28) 

The dynamics of H(t, Xi,X2) is obtained by using Ito's formula (see Appendix [D|) to differ- 
entiate Eq. (1271) with respect to t and then eliminating ^ with the help of Eq. The 
result is 
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;H{t, xi, X2) =Ds + ^ ) H{t, xi, X2) - DgH{t, xi, X2)6{xi - X2) - 2(/ii2 + ^2i)H{t, xi, ^2) + 



dt ' ' \dx\ dxl 



2/i21 + (/il2 - fi2l)[F{t,Xi) + F{t,X2)] + -[H{t, Xi,Xi) + H{t,X2,X2)]- 
s{2f{t,X,)[l - f{t,Xi)]f{t,X2)) - s{2f{t,X2)[l - f{t,X2)]f{t,Xi)). 

(29) 



Equations (!28l) and (|29l) agree with the ones derived in Ref. 25|] in the hmit of no mutations 
considered there. 

From Eq. fl29|) . one can see that the hierarchy of the moment equations does not close 
unless selection is absent. Similar to the well-mixed case, the correlation functions for 
neutral models with and without mutations can be found analytically, see Sees. |V] and HV] 
but different methods are required to analyze the dynamics in the presence of selection, 
see Sec. I VII To simplify the analysis, we consider well-mixed, spatially homogeneous initial 
conditions. Then F is only a function of t, and if is a function of t and x = Xi — X2- 
With these simplifying assumptions, the equations of motion for F{t) and H{t,x) take the 
following form: 

^ = /^2i - (/ii2 + fi2i)F{t) + '-H{t, 0), (30) 



-H{t, x) =2Ds-^H{t, x) - DgH{t, 0)6{x) - 2(/ii2 + fi2i)H{t, x) + 

2/221 + 2(/ii2 - fi2i)F{t, x) + sH{t, 0) - 2s(2/(t, xi) [1 - f{t, xi)]f{t, X2)). 

IV. NEUTRAL MODEL WITHOUT MUTATION 

We start the analysis of the one-dimensional stepping stone model by considering neutral 
alleles that do not mutate. In practice, this means N'^fli2, N'^fl2i ^ 1 and N'^s <^ 1 (as 
we show below). Although these assumptions are not always realistic, they help to clarify 
the role of genetic drift in a spatial context. In addition, neglecting mutations is a good 
approximation on time scales shorter than the waiting times for the mutations jji2 ^^(^ l^2i- 
Under these assumptions, F does not change, F{t) = Fq, and Eq. fl3T|) reads 
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-H{t, x) = 2D,—H{t, x) - DgH{t, 0)6{x). (32) 

Equation (l32l) can also be derived by tracing the lineages of organisms backward in time. 
The spatial heterozygosity, H{t,x), is the average probability of sampling two different in- 
dividuals chosen at time t from demes separated by distance x. As we trace the hneages 
of the two sampled organisms backward in time, the lineages diffuse in space due to mi- 
gration and when they are at the same point they have a chance to coalesce, in which case 
the sampled organisms must be identical because they have a common ancestor. Such a 
coalescence event changes the probability of being different from H to 0, and acts like a sink 
at X = 0. The first term on the right hand side of Eq. ( l32l) describes the diffusion, and the 
second term describes the coalescence. Since this argument is valid for an arbitrary number 
of alleles, Eq. (!32l) is valid for an arbitrary number of spatially diffusing neutral alleles. See 
Appendix [E] for a more detailed discussion of the g-allele problem, with q > 3. 

To better understand the microbiology experiments on neutral alleles by Hallatschek et 
al. 5|, we consider uncorrelated initial conditions -F(O) = Fq and if(0,x) = i^o, where Fq is 
the fraction of allele one and Hq = 2Fq{1 — Fq), which is the heterozygosity of a well-mixed 
population with the frequency of allele one equal to Fq. For these initial conditions, Eq. (!32|) 
is solved in Appendix \M The results are 



Hit, .).H,- I it'^=== erfc | e-. (33) 




i/(i,0) = i/oerfc . ens, (34) 



where eiic{y) is the complementary error function. 

The spatial heterozygosity at vanishing separation, H{t,0), is particularly interesting 
because it indicates the degree of spatial segregation: if H{t,0) <ti 1, then, locally, the 
demes are fixed to one of the two alleles. From Eq. (!34l) . one can see that, for t ^ SDs/D^, 

^^(i,0)=/Jo(^^j +0{t-'/'), (35) 

which means that at long times one of the alleles reaches fixation locally. Therefore, we 
see that the spatial model we are considering is consistent with the experiments by Hal- 
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FIG. 4: (Color online) Solutions of Eq. (I32p at various times, given H{0, x) = ^. Time and distance 
are measured in untis such that Ds = 1 and Dg = 1. Time increases from the top curves to the 
bottom curves. Note the statistical reflection symmetry, H{t, x) = H{t, —x) 

latschek et al. jsl (see Fig. [1]) because it predicts the formation of domains (the regions of 
local fixation). Thus, similar to the well-mixed model considered in Sec. [Tll one of the alleles 
reaches fixation locally with fixation time rj = 8Ds/(7rD|) ~ A^^. But not only is this 
fixation time proportional to A^^, instead of A^, the functional form of heterozygosity decay 
different: instead of a rapid exponential decay, the spatial model shows a slow algebraic 
decay of local heterozygosity. These results agree with the previous works on population 
genetics by Malecot [13] and Nagylaki [24]. Local fixation and t"^^"^ decay of local heterozy- 
gosity have also been found in the voter model [12] , which corresponds to the stepping stone 
model with = 1. 

The characteristic demixing time can also be estimated by the following scaling argument: 
the effective population size at time t in the coarsening process is Neff{t) ~ nQy/tD^, where 
the population density tiq ~ N/ a. Upon recalling that the fixation time in zero dimensions 
is Tf = NeffTg, and solving self-consistently for tj, we have tj ~ DgN^Tg/a^ ~ Ds/Dg ~ 
N\. 

Another important characteristic of H{t, x) is the length scale over which H{t, x) changes 
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from its minimum to its maximum values. Figure H] plots Eq. (l33ll and shows the spatial 
variation of H(t, x) at different times. One can see that the spatial heterozygosity is reduced 
near the origin due to the local fixation, but H{t,x) rises to its initial value Hq at large x, 
where the alleles remain uncorrelated. After the domains form, i.e. for t ^ SDs/D"^, 
this change from H{t, 0) to H{0, x) happens on a length scale that is set by the average 



size of the domains £, which is proportional to the diffusion length ^j2Dst^ as follows from 
Eq. fl33|l . Since this characteristic length scale changes with time, it is convenient to rescale 



distances: x = xj ^2Dst. Upon using Eq. (135|) to simplify Eq. (133!) . we see that H{t,x) 
approaches a nontrivial limit in terms of x as time goes to infinity: 

,3a) 

which agrees with the known results for the voter model 

A more precise evaluation of the domain density and hence an average domain size i{t) 
can be obtained by differentiating H{t,x), as shown in Appendix [Bl From Eq. flB4l) . we 
know that £{t) = jj^^^, so using Eq. (135|) we see that 



2/o(1f4o)' 



which is consistent with the analysis of Ref. [2 



Note that the genetic diffusion con- 



stant Dg ~ drops out. As discussed in Ref. [26|, at large times the only dynamics left is 
the motion of the domain boundaries; with neutral alleles, these walls behave as annihilating 
random walks, and the average domain size can be easily calculated. 

Equations ( IMl) and ( 1371) suggest that the processes driven by the genetic drift slow down 
with time because the logarithmic time derivatives of H{t, 0) and i tend to zero as time goes 



to infinity. In the annihilating random walk picture of Ref. [26l], annihilations become rarer 
and rarer as the coarsening progresses. A more direct measure of genetic drift, which is also 
interesting from the biological point of view, is the fluctuations of the total fraction of, say, 
the first allele f{t) in a finite population of length L. We define f(t) as 

1 '•^ 



fit) = l I f{t,x)dx, (38) 







and compute its variance z/(t) to characterize its fluctuations. 
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Upon integrating Eq. (l2^ over x with s = /ii2 = fi2i = 0, we obtain the equation of 
motion for f: 



Jt = j£ ^jDJ{t,x)[l-f{t,x)nt, x)dx, (39) 

where the spatial diffusion term vanishes after integration by parts provided periodic or 
Newman boundary conditions are imposed. Upon noting that (f) = F = const (the Ito 
interpretation of the noise T(t,x) is crucial here) and defining 

^ = if) - (f)^ (40) 

one finds immediately that ^ = ^(f^)- To evaluate the time derivative, we use the rules of 
the Ito calculus, sketched in Appendix [D| and find 



^ = ^J j ^jDJ{t,xr)[l - f{t,Xr)]^DJ{t,X2)[l - /(t,X2)](5(xi - X2)dx^dx2, 

(41) 

where the delta function comes from averaging over the noise and using Eq. (l25l) . From 
Eq. gH it follows that 

u{t) = ^I^H{t\0)dt\ (42) 

where we assume z/(0) = 0. Hence, we know z/(t) exactly because H{t, 0) is given by Eq. (IMI) . 

For small times, t ^ SDg/Dg, the variance grows linearly with time. For large times, we 
can use the asymptotic expansion of H(t,0) given by Eq. (135!) to calculate i^{t). The result 
is 

.W.^^ofi^^V (43) 



J^L \DgL 



Equation (H3|l is consistent with Ref. [16|, and we immediately conclude that the standard 
deviation A(t) = \/i4t) grows as t^/^ for large times! This important result for g- alleles 
is re-derived in Appendix [E] by approximating the dynamics of the domain boundaries by 
annihilating random walks. Thus, f performs a subdiffusive random walk, and the action of 
genetic drift on the global frequency f{t) becomes weaker with time. Equation (H3l) is valid 
only for t <^ ^ because it relies on Eq. (135!) . which is valid for an infinite population, and 
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should break down at times that are long enough for a domain boundary to diffuse from one 
end of the population to the other. Using Equations ([8]) and fj42|) . one can also calculate 
the behavior of the global heterozygosity Ti.{t) = L^^ H(t,x)dx, i.e. the probability to 
sample two different alleles from the population regardless of their spatial locations: 

2F„(l- F„) - = 3. £ H(t', 0)dt' = + O (^) . (44) 



where the second equality requires t <^ for the reasons mentioned above. In the opposite 
limit t ^ the global heterozygosity T-C{t,x) obeys zero-dimensional dynamics of a well- 



mixed population with an effective population size L/{2Dg) as shown in Ref. 24l |. 

The local heterozygosity and average domain size can be obtained from experiments on 
microbial spreading like the one shown in Fig. [H If the data are sufficiently precise, Eqs. f p5|) 
and ( 1371) could be used to extract Ds and Dg from the experiments. Since Dg ~ 
extracting Dg from experimental data determines the effective deme size for the equivalent 
stepping stone model. Dg can be obtained from the diffusion of individual domain boundaries 
or z^(t). These two parameters completely determine the neutral dynamics without mutation 
and play an important role when selection or mutation is present. 

V. NEUTRAL MODEL WITH MUTATION 

While on short time scales mutation can be neglected, it is the long time scales and the 
patterns of genetic diversity created by mutations that are of particular interest in population 
genetics. Noticeable mutations also arise in microbiology experiments like those in Fig. (H 
especially if mutation rates are enhanced by DNA damaging chemicals or radiation. In this 
section, we extend the results of Sec. |IV]by allowing for nonzero mutation rates between the 
two alleles. We assume as before statistically homogeneous initial conditions, and note that 
the dynamics of the dynamics of the one and two-point correlation functions is then given 
by 

= 1^21 - (1^12 + l^2i)F{t), (45) 
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-H{t, x) =2D,—H{t, x) - DgH{t, x)6{x) - 2(/ii2 + fi2i)H{t, x)+ ^^^^ 

where F{t) = {f{t,x)) is independent of x. The equation of motion for F in the spatial 
model is exactly the same as Eq. ( ITTl) . which describes the well-mixed-population model. 
Therefore, F relaxes to its equilibrium value, F{oo) = ^7^^^^^? [see Eq. f|T3l) ] exponentially 
fast with time constant (/ii2 + f^2i)^^ ^ t^- The similarity to the nonspatial model is not 
surprising because neutral mutations occur equally likely at any point within the population, 
regardless of its spatial structure. The dynamics of H{t, x) is, however, more complicated 
because both mutation and genetic drift determine the behavior of the spatial heterozygosity. 
The stationary solution of Eq. fl46l) reads 



(_ / M12+M21 I J. 
1 _ " ^ I (47) 

Equation (H7|) agrees with the solution by Kimura and Weiss [6], which was obtained in the 



discrete space and time limit. One can see that, for x ^ y /Ji2+At2i ' spatial heterozy- 
gosity approaches 2F(oo)[l — F(oo)]. Thus, mutations cause the frequencies of allele one to 
eventually become uncorrelated at large separations. At shorter distances, however, there 
are correlations, and H{oo,x) < H{oo, 00) for all x < 00. Note, in particular, that 

2£M11-^^ (48) 



4 V Ds(/il2+/i2l) 

Note also that, for small mutation rates, the heterozygosity is proportional to /ii2 + fi2i [see 
Eq. f[T^ ] in a well-mixed population, but the local heterozygosity in a one-dimensional 
population is proportional to VyUi2 + /^2i whenever Tg{fii2 + fi2i) / (N'^DsTg), which is 
a consequence of weaker genetic drift in one dimension ^. 

When if (00, 0) ^ H{oo, 00), the population is segregated into domains of different allelic 
types. Upon invoking Eq. (1B4I) we obtain the following average domain size: 



^ In population genetics, population structure and spatial correlations are often reported via Fgt 



H{oo^ 0)/iJ(cx), 00), which can readily be obtained from Eq. 
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DgfIl2fJ'21 \ 4:\l Dsiflu + fJ'2l) J 2/ii2Ai21 

in the limit discussed above. This result, together with Eq. (fT3|) . can be used to extract the 
mutation rates from experimental data. 

We can also determine how fast H{t, x) reaches its stationary value. Since the het- 
erozygosity cannot be in equilibrium unless the frequency of the alleles has equilibrated, we 
assume for simplicity that F{0) equals its stationary value. Then, the deviation of the spa- 
tial heterozygosity from its long time equilibrium value H{t,x) = H(t,x) — H{oo,x) obeys 
the following equation: 

-H = '^Ds-^H - D,H5{x) - 2(/ii2 + ^^2l)H. (50) 

Equation ( 150|1 can be further simplified by the change of variables H = e~'^^^^^^^^'^^^H , 
which leads to 

Wt^ = '^^^^2H-DM^)- (51) 

Since Eq. f lSTl) is identical to Eq. fl32l) . we conclude that, at long times, the difference 
between H{t,x) and the stationary solution decays as C^^^^^e'^^'^i^+^^i)*^ where C is a 
constant. Thus, apart from an algebraic prefactor (and a nontrivial spatial dependence), 
the dynamics of H{t, x) is essentially the same as in the well-mixed case. 

In this section, we considered a model with only two alleles; however, in many circum- 
stances, an infinite alleles model is more appropriate. The infinite alleles model is briefly 
discussed in Appendix O Some results forg-alleles, 2 < g < oo, are discussed in Appendix |El 

VI. SELECTION 

Unlike the neutral models with spatial diffusion and mutation discussed above, the one- 
dimensional stepping stone model with selection cannot be treated analytically because the 
hierarchy of moment equations does not close. We briefly examined three closure schemes: 

(2/(t,xi)[l - /(t,xi)]/(t,X2)) ^ Hit,0)F{t), (52) 
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{2 fit, xi) [1 -fit, x^)]fit, X2)) ^ 2F(t) [1 - Fit)] [1 - 2F(t)] - [1 - 2F(t)]/J(t, x) +Hit, O)F(t), 

(53) 

and 



(64) 

The first scheme is a simple factorizatio n ap proximation: the second scheme, which assumes 



small fluctuations, is due to Nagylaki |25|]; and the third scheme, which provides a good 
approximation for some diffusion limited reactions, was proposed in Ref . . Unfortunately, 
none of the schemes describe the behavior of the system correctly. The progress can be made, 
however, for some initial conditions in two limiting cases of strong selection ^Mf- ^ 1 and 

weak selection <^ 1. For the rest of this section, we include spatial diffusion and genetic 

3 

drift, but neglect mutations, which is justified on short time scales. 

First, let us consider the initial condition /(O, x) = 1 — 6ix), where 6'(x) is the Heaviside 
step function. This initial condition specifies just one domain boundary, which, for any 
positive s, undergoes Brownian motion with a drift to the right. This is a good description 
of an expansion of a new advantageous mutant spreading through the population. In the 
strong selection limit (s ^ D^/Ds), R. A. Fisher [9| found that the sharp boundary above 
broadens to a width of order ■>/ Dg/ s, and the velocity of the genetic wave is given by 



Vs = 2V^. (55) 
When, in contrast, selection is weak compared to genetic drift, Doering, Mueller, and 



Smereka 



recently found that the velocity is given by, 

2sD, 



(56) 



When the population contains multiple domains, the domain walls bordering a favorable 
genetic variant ( "allele one" ) expand to engulf the regions occupied by the more deleterious 
allele. 

Another interesting initial condition is /(0,x) = Fq = const., i.e. the population is 
initially well-mixed. This scenario, for example, describes the quasi-one-dimensional strip of 
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pioneers advancing at the front of a two-dimensional population wave that originated from 
a well-mixed ancestral population and is propagating in the region where one of the alleles 



has higher fitness (see Ref. [26|). If selection is strong enough, then allele one ("the preferred 
variant") takes over the population before spatial correlations have time to appear. To see 
this, note from Eq. (1181) that allele one wins locally on the time scales of s~^, but, from 
Eq. the time for spatial correlations to appear is on the order of Dg/Dg, which is much 
larger than in the strong selection limit. Thus, the behaviors of one-dimensional and 
well-mixed populations are similar when s 3> D^/Dg = a? / {r'^N'^Ds). 

In the limit of weak selection, however, spatial correlations appear before allele two is 
eliminated, which changes the dynamics of the system. Qualitatively, we can divide the 
selective sweep into two stages. During the first stage, the effects of selection are negligible 
and spatial segregation occurs as described in Sec. [IVl During the second stage, the domains 
of allele two shrink at each end with wall velocity given by Eq. (!56|) . and the stochastic 
motion of the domain boundaries can be neglected. The crossover time between the stages 
occurs when the diffusive displacement of the walls is of the same order as their deterministic 



displacement, that is when Vwt = ylDgt- Thus, the crossover time r* is on the order 
of Dg/vlj, which can be expressed as Dg/{s'^Ds) with the help of Eq. fl56|) . Then, from 
Eq. (!37|) . the average domain size at the crossover i* is on the order of Dg/{sHQ). 

The dynamics during the second stage depends on the probability distribution of domains 
of size rj, Pd{f]) at time r*. For annihilating random walkers, which are a^ood approximation 
to domain boundaries during the first stage, Bramson and Griffeath [29] proved that Pdij]) 
has exponential tail for large 77 of the form e"''"'/^*, where 7 is a number of the order 1. Since 
each domain shrinks with velocity 2f^, the fraction of allele two can be expressed as 



7(i7 + 2nmt) 



1 — F{t) oc / rje dr] oc e , (57) 

Jo 

where A is a number of order 1. From Eq. (!57|) it follows that, as in the well-mixed case, the 
selective sweep is exponentially fast, but the time constant of this process is proportional 
to rather than s~^. 

The analysis leading to Eq. ( 1571) can be generalized to an arbitrary initial probability dis- 
tribution, Pdivi), provided the dynamics is dominated by selection and genetic drift. For ex- 
ample, if a population initially in equilibrium with respect to mutations and genetic drift (see 
Sec. |V]) is affected by an abrupt environmental change that makes allele one advantageous, 
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then the shift to the new equihbrium occurs exponentially fast with a time constant propor- 

1/2 

tional to ^^^^ 2+^1^21)^^^ ' assuming Pdiv) has exponential tail of the form e"'^'''/^, where 7' is 
a constant of order unity, and i is given by Eq. fl49p . 

We make two important observations based on the results of this section. First, the 
temporal dynamics of the one-dimensional stepping stone model with selection can depend 
strongly on the initial conditions. Second, the results in the weak selection limit are related 
to the results in the strong selection limit by the transformation s s'^Ds/D'^ at least up to 
a numerical factor. The second observation suggests that, while data can be naively fitted 
to a well-mixed-population model, the fit in fact gives the "renormalized" value of s instead 
of the "bare" one. 

VII. SIMULATIONS 

In Sees, mil HVl IVl and |VT] we reviewed and extended the theoretical analysis of the 
one-dimensional stepping stone model. This model, while of great theoretical interest, relies 
on a restrictive set of assumptions including large deme sizes and slow diffusive migration. 
The recent experiments by Hallatschek et al. [5], on the other hand, were carried out with 
bacterial fronts that were only a monolayer thick; therefore, demes consisted of only a few 
microbes. Moreover, depending on the microorganism, nearby demes can exchange a signif- 
icant fraction of cells in each generation. In this section, we discuss numerical simulations 
not subject to these restrictions and compare them with the theoretical predictions. 

We simulate L organisms arranged on a line and labeled by an integer /, / = 1,2,...,L. 
Each organism can be either of allelic type 1 or allelic type 2. During even generations, 
the offspring at site / comes from an organism at either site / — 1 or site /, whereas, during 
odd generations it comes from either site Z or Z + 1. The simulations embody the process 
illustrated in Fig. [2], laid out on a triangular lattice in space and time. Periodic boundary 
conditions are imposed at the left and right ends of the population. Let 12 ^ 2 refer to 
the event that the offspring has allelic type 2, while one of its possible parents has allelic 
type 1 and the other has allelic type 2, etc. The transition probabilities, which depend on 
the states of the possible ancestors, are then given by 
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11^2: /2i2, 

22 ^ 1 : /i2i, 

22 ^ 2 : 1 - /i2i, 

12^1: (1 -/iia) + ^21, ^ ^ 

1 + 5 1 — 5 / 
12^2: __^i2 + ^(l-/i2i), 

1 + 5,. 1 — 5_ 

21^1: (1 -/ii2) + /i2i, 

1 + 5 1 — 5 , 
21^2: __/,,2 + ^(l-^2i). 

The event 12 — 2 can happen either if allele one was selected for reproduction [probabil- 
ity (1 + 5)/2] and it mutated (probability /ii2) or if allele two was selected for reproduc- 
tion [probability (1 — 5)/2] and it did not mutate (probability 1 — /i2i). Other transition 
probabilities are obtained analogously. Thus, the system we simulate is very similar to the 
voter model 1^, equivalent to population genetics models with = 1. However, to make 
the calculation faster, we use discrete generations rather than exponentially distributed wait- 
ing times until reproduction. We found no significant differences in dynamics between the 
voter model and the model used here. 

First, we simulate the neutral model without mutations. To illustrate the similarities and 
differences between the one-dimensional stepping stone model and a growing population 
with an undulating front, we also simulate a flat population wave in a two-dimensional 
habitat;^ both are displayed in Fig. [21 Our model with an undulating front is the same as 
in Ref. [11], but we use a triangular grid instead of a cubic one. Figures [5] and [6] show how 
the average number of domain boundaries decrease with time; the insets show the mean 
square displacements of the boundaries as a function of time. The simulations confirm that 
the domain boundaries perform a diffusive random walk for the (1 + l)-dimensional model 
with a fiat front, but they move superdiffusively for the undulating front (Eden) model, in 



agreement with Ref. ll|. In this case, the number of domain boundaries decays faster for 
the undulating-front model compared to the (1 + l)-dimensional fiat front model that, as we 
show below, most closely tracks the prediction of the one-dimensional stepping stone model. 
A single run of a simulation is shown in Fig. [TJ and the spatial heterozygosity averaged 
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FIG. 5: (Color online) The number of monoallelic domain boundaries as a function of time in 
the undulating front model. The simulation of 100 demes averaged over 100 runs is plotted in 
blue (dots), and the theoretically expected decay of the number of boundaries as (see Ref. 11|) 



is plotted in red (solid line). Inset: Log- log plot of the mean square displacement of monoallelic 
domain boundary end-points as a function of time in the same set of simulations as the main plot. 
The blue dots are the simulation data, and the solid red line is the expected slope according to 
Ref. ill|. 

many runs is shown in Fig. [HI Note that coarse-graining is required to properly represent the 
time evolution of H{t, 0). Figure [3 shows that H{t, x) for large t can be described well by the 
limiting shape of spatial heterozygosity given by Eq. (!36|) . To properly represent H(t,0), 
we artificially define demes of a larger size by grouping M neighboring individuals into 
one deme. From a theoretical point of view, this procedure is similar to the formation of 
Kadanoff block spins as in renormalization group methods [s^ whereas, from the point of 
view of population genetics, this procedure is similar to the methods of collecting data from 
a dispersed natural population. In field studies, scientists do not typically sample every 
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FIG. 6: (Color online) The number of monoallelic domain boundaries as a function of time in the 



flat front model. The simulation of 100 demes averaged over 100 runs is plotted iri 
and the theoretically expected decay of the number of boundaries as (see Ref. 



jlue (dots), 



26|]) is plotted 



in red (solid line). Inset: Log- log plot of the mean square displacement of monoallelic domain 
boundary end-points as a function of time in the same set of simulations as the main plot. The blue 
dots are the simulation data, and the solid red line is the expected slope according to Ref. l3, [3] 
and Eq. (fHTll . 

single individual; instead, they often divide the habitat into patches and sample a repre- 
sentative number of individuals from those patches. To summarize, we keep the dynamics 
of the simulation exactly the same, but define the spatial heterozygosity on the demes of 
size M rather than on single individuals. We found that the local heterozygosity has the 
form H{t,0) = (3{M)t-^/^, as predicted by our analysis of the stepping stone model, for 
all M studied (1 < M < 64). From Eq. (l35l) . we expect that [3 oc M~^; this expectation 
was also confirmed by our simulations. 

As discussed in Sec. llVt the total fraction of allele one f(t) fluctuates in an unusual way 
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FIG. 7: (Color online) Simulation of microorganisms growing along a cylinder, (a) Cross sections 
of the cylinder at the inoculation time (bottom) and at two later times. The two genotypes are 
shown in red and green. Initially, both alleles are present in almost equal proportions, and the 
two later cross sections show the results of the coarsening process, (b) The spatially averaged 
heterozygosity [defined in the sense of Eq. (|Bip but without taking the limit L — > oo] for the three 
cross sections shown in a; the separation / is the shortest distance between two points around 
the cylinder. The colors correspond to the colors used is a. At inoculation, the heterozygosity 
fluctuates around 1/2 since the two phenotypes have equal probabilities of occupying any site. The 
only exception is the site / = 0, where the heterozygosity is zero automatically because we only 
allow a single microorganism per site. After 1024 generations short range correlations are clearly 
visible and after 4096 generations one can easily relate the features of the heterozygosity function 
to the sizes of the sectors in the corresponding cross section. The wiggles are eliminated when 
averaged over many realizations, as shown in Fig. [8l 

with time. Figure [TUk shows examples of these remarkable variable-step-length random 
walks. The fluctuations of f(t) obey Eq. and grow sub diffusively, as shown in Fig. [TUb . 
We also find good agreement between the theory and the simulations in the presence of 
mutation for all values of M studied. Figure [11] shows the stationary spatial heterozygosity 
for M = 1. Finally, we studied selective sweeps in an initially well- mixed population. Fig- 
urefT^ confirms the prediction from Eq. flFTI) that F oc (1 — e~"*), and Fig. [T2b confirms the 
result of Sec. |Vl]that, for strong genetic drift, the effective extinction rate a is proportional 
to s^. 

Numerical results for three neutral alleles are presented in Appendix [El 
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FIG. 8: (Color online) Effects of coarse-graining of the time evolution of the heterozygosity H(t, I). 
The curves are averages of 100 model runs, (a) The (1 + l)-dimensional model with a flat front, 
where each deme hosts only one organism. Consequently, the heterozygosity at / = is zero at 
all times. Note the qualitative similarity to the analytical results at larger distances except for 
boundary effects, (b) Similar runs for the model except that the organisms have now been grouped 
into demes of size 5 before calculating the heterozygosity. This coarse-graining gives qualitative 
agreement with the analytical results even at / = 0. 

VIII. CONCLUSIONS 



Fluctuations due to sampling error during reproduction significantly affect the evolu- 
tionary dynamics of quasi one-dimensional populations, e.g. two-dimensional populations 
undergoing range expansions. These fluctuations lead to the genetic demixing illustrated in 
Fig. [131 where an initially well-mixed population of alleles "phase separates" into monoal- 
lelic domains. The transition is somewhat analogous to spinodal decomposition in physics 
and material science 3l|, but is also markedly different. In particular, unlike conventional 



demixing phase transitions in statistical mechanics, genetic demixing occurs only in a low 



number of spatial dimensions d {d < 2) 



21 



31| . The dependence of genetic demixing on 



the number of spatial dimensions d is illustrated by the decay of local heterozygosity in the 
absence of selection and mutation. For long times, the functional form of the decay is given 

by y. 
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rescaled distance, a; 



FIG. 9: Comparison between the analytical prediction for the limiting shape of H{oo,x) and 
the simulations with a flat population front. The continuous curve is formed by the data points 
representing H(t, x) for several times between t = 2 ■ 10^ and t = 4 • 10^ plotted in the rescaled 
coordinates x, and the circles represent the theoretical prediction of the limiting shape of the spatial 
heterozygosity, see Eq. (j36p . The data are obtained in a simulation of 3200 individuals for 4 • 10^ 
generations with averaging over 500 realizations. At t = 0, each site was assigned either allele one 
or allele two with equal probability. 



-t/i2NTg) 



H{t,0) ~ < 



d = 0, 
d=l, 
d = 2, 
d>2. 



(59) 



l/ln(t) 
const 

Note that = 2 is the critical dimension. 

Here, we have shown that the one-dimensional stepping stone model has very different 
dynamics compared to the standard well-mixed-population models used in genetics. Most of 
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FIG. 10: (Color online) Genetic drift in a finite population, (a) The total fraction of allele one f(i) 
versus time in four single runs of the neutral model with a flat front. Here L = 1000, and there 
are no mutations, (b) The standard deviation of the frequency of allele one A(f), shown with blue 
dots, is obtained from 220 realizations of the simulations described in a. The red solid line shows 
the slope expected from. Eq. (jl3]) . 

the differences arise because, in the spatial model, populations segregate into monoallelic do- 
mains. As a result, genetic drift and selection can only act at the boundaries of the domains, 
which slows down the dynamics of the model. In particular, we found that, in the neutral 
model without mutation, fixation occurs exponentially fast in a well-mixed population, but 
the decay of heterozygosity is algebraic in the spatial model. Genetic drift in the population 
as a whole becomes weaker with time, as spatial diffusion causes the effective population 
size to increase. We also found that the standard deviation of the total fraction of one of the 
alleles (in the absence of selection and mutation) increases sub diffusively as t^^^. Selective 
sweeps also occur more slowly in the spatial model: for weak selection, s <^ D^/Ds, we 
found that the time constant of the selective sweep is quadratic in s in the spatial model, 
but it is only linear in s in the well-mixed-population model. The effects of mutation do not 
differ as dramatically in spatial and nonspatial models, but spatial stepping stone model 
reveals nontrivial spatial correlations and predicts a different value for the local steady state 
heterozygosity, proportional to y/^u + /i2i for small mutation rates, compared to /xi2 + fJ.21 
in the well-mixed-population model. 

Our main conclusion is that the data from natural populations may not always conform 
to the predictions of the well-mixed-population model and, even when it does, the estimated 
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FIG. 11: Equilibrium between mutation and genetic drift in the absence of selection. Comparison 
between the analytical prediction for the steady state heterozygosity H{oo,x) and the simulations 
with /ii2 = 10~^ and fl2i = Ai2- The dots represent the results of the simulation, and the circles 
represent the best fit of theoretical result given by Eq. (f47|l to the data. The data are obtained in a 
simulation of 3200 individuals for 2- 10^ generations with averaging over 100 realizations. At t = 0, 
each site was assigned either allele one or allele two with equal probability. 



parameters from the model may not be biologically meaningful. The spatial model contains 
an important additional parameter, the spatial diffusion constant parameter Ds, which enters 
into many of the predictions. For example, the timescale for local fixation is given by Dg/D"^ 



rather than NTg [see Eq. fl35l) ] and, for weak selective advantage, s is replaced by s'^Dg/Dg, 
see Sec. IVII Moreover, as we saw in Sec. I VII I the timescale of fixation depends on the 
partitioning of the population by the experimenter into measurement sites. Thus, care must 
be taken when interpreting the data from the natural populations. Finally, well-mixed- 
population models and experiments without spatial resolution do not account for spatial 
correlations, which contain important information about the population. 



35 



400 600 800 1000 
t 

I -3 

ln(s) 

FIG. 12: The effective extinction rate a versus s in tlie limit of weak selection. The dots are the 
data from the simulation, and the line has the slope equal to 2, that is the expected slope from 
Eq. (|57|) . The data supports a oc s^.The values of a are obtained from graphs like the one shown 
in the inset. Inset: ln(l — F) versus t for s = 0.06. The circles are the actual data points, and the 
line is the best least squares linear fit. The simulation confirms exponentially fast fixation. The 
data are obtained in a simulation of 1600 individuals for 6000 generations with averaging over 100 
realizations. At t = 0, each site was assigned either allele one or allele two with equal probability. 
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FIG. 13: (Color online) Illustration of spinodal-decomposition-like genetic demixing in a one- 
dimensional population, (a) Initially well-mixed population with red and green colors labeling 
different genotypes, (b) The same population several generations later. The frequency of one of 
the alleles is now oscillating between and 1 because the population segregates into monoallelic 
domains. 
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APPENDIX A: SOLUTION OF THE NEUTRAL MODEL WITHOUT MUTA- 
TION 



In this appendix, we solve Eq. (15^ subject to initial condition H{0,x) = Hq. It is 
advantageous to first solve a simpler equation: 

§-tH = ^^^-^2H-bmx), (Al) 

where b{t) is an arbitrary function of time. Equation ( lAip is a standard diffusion equation 
with a sink term, and it can be readily solved in the Fourier domain. The result is 

/•* b(t')e »Os{t-t') 
H(t,x)=Ho- / dt'-^^^ (A2) 

Note the convolution of b{t') with the diffusion propagator. Now, we impose a self- 
consistency condition b{t) = DgH{t,0), which leads to 

This is Abel's integral equation of the second kind, canonically written as 

y{x) + X r^^ = g{x), (A4) 

J a -t 



where g{x) is a known function. The general solution of Eq. (IA4p is given in Ref. 33], and 
it reads 

y{x) = G{x) + ttA^ J e"^'(^-*)G'(t)rft, where 

r^r ^ r ^ ^ ["'^ 9{t)dt ^^^^ 

J a -t 

Equations ([33D and ^ follow from Eqs. (|A2|, (jAS]), 0, and (lASD . 



APPENDIX B: AVERAGE DOMAIN DENSITY FROM THE SPATIAL HET- 
EROZYGOSITY H{t,x) 

In this appendix we derive the relationship between the spatial heterozygosity, H(t,x), 
and the average domain density nd{t). From n^, we can get a domain size by defining i = n^^. 
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The result for the domain density is vahd for an arbitrary number of alleles, so in this 
appendix we use a broader definition of H{x,t) as the average probability of sampling at 
time t two different alleles from two demes a distance x apart. We assume that the domains 
have formed, and they are on average much larger than the boundary regions. 

Let h{t, xi,X2) equal to 1 if both xi and X2 are occupied by organisms in different allelic 
state and otherwise. To compute i, we use an alternative definition of H{t, x) with ensemble 
average replaced by space average: 



H{t, x) = hm J [ hit, e, e + x)di, (Bl) 

where we assume periodic boundary conditions. Let us compute H(t,x + 6x) — H(t,x) = 
lim^^^oo X /o^[^(^'^'^ + ~ ^(^) 'C) Ol'^'C foi' ^x small compared to typical domain size but 
large compared to the deme spacing a. Now expand both sides in 6x. At the lowest order 
in Sx, each domain boundary contributes 5x to the right hand side; therefore, -^H{t, +0) 
equals the density of the domain boundaries. Upon defining the average domain size i{t) as 
the inverse of the domain boundary density, we obtain the following relationship: 

This relation is analogous to the one derived in Ref. 33|. 
We can further simplify Eq. flB2l) by observing that 

dx = ' ^^^^ 

which follows from integrating Eq. ( l3Ti) or Eq. ( ICll) with respect to x from — e to e, < e <^ 1, 
and noticing that H{t, x) is an even function of x. The final result then reads 

It should be emphasized that this result is only valid in the limit of very large domain sizes 
compared to the boundary regions, which means H{t,0) <^ 1. Therefore, the leading term 
in H{t, 0) is sufficient at this level of approximation. Note that Eqs. IB2I and IB4I are valid 
in the most general settings considered in this paper, i.e. in the presence of genetic drift, 
migration, selection, and mutation. 
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APPENDIX C: INFINITE ALLELES MODEL 



In this appendix, we extend the analysis of the stepping stone model with mutations 
presented in Sec. |V] to the infinite alleles model. The infinite alleles model assumes that 
every new mutation creates a new allele, which is a good approximation for genes encoded 
by a large number of nucleotides because the number of all possible mutations is much larger 
than the number of all possible back mutations [l7]. The equation of motion for H{t,x), 
which we interpret as the average probability of sampling two different alleles from demes 
X apart, can be derived by following two lineages backward in time, as done in Sec. IIVI In 
the presence of mutation, the right hand side of Eq. (152]) should contain an additional term 
describing the rate of increase of H{t, x) due to mutations in both of the lineages. Because, 
in the infinite alleles model, a mutation changes the probability that the organisms have 
different alleles from H to 1, that is by 1 — H, the new term is 2/i(l — H), where /i is the 
mutation rate that is assumed to be the same for all types of mutations. Thus, Eq. fl32|) 
becomes 

-H = '^D^—H + 2^(1 -H)- D,H6{x) (CI) 

for the infinite alleles model (compare Eq. H6l) . The stationary solution of Eq. fICll) is given 

by 

g(oo,x) = l- ^ (C2) 
1+1,/^ 



At large separations, H{oo, x) approaches one, which is consistent with the infinite number of 
alleles. Locally, H{oo, 0) = (1 + \\J '^^)~^: and if H{oo, 0) ^ 1 the population is segregated 
into domains containing only one allelic type. The average size of such domains follows from 
Eq. (IB4]): 

1 / CM ID. 

/i 

where the last equality follows from the assumption that H{oo,0) ^ 1. The approach to 
the stationary state can be either obtained by methods of Appendix |A] or by the change of 
variables H{t,x) = H{oo,t) + e-'^^'*H{t,x), which reduces Eq. ([Cl]) to Eq. ([32]). The result 
is that the slowest decaying mode vanishes as Cw-^e"^'^*, where C is a constant. 
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The infinite allele model and Eq. (1C1|) has been analyzed before in Refs. 13|, |2J], where 
the stationary solution and the long time approach to the equilibrium were calculated. Our 
results are consistent with their findings. 



APPENDIX D: THE ITO CALCULUS 

In this appendix, we briefly discuss the Ito calculus. Our presentation relies on Refs. [Ji 
and [22^, which can be consulted for a more extensive presentation. For simplicity, we 
only consider nonspatial stochasic differential equations, but the results can be extended to 
spatial problems straightforwardly. 

Let us analyze the following stochastic differential equation, which includes Eq. (!20l) as 
a special case, 

^=u{ij)+g{^)T{t), (Dl) 

where r{t) satisfies Eq. ( |2T|) . and ujlip) and g{ip) are arbitrary continuously differentiable 
functions. From the point of view of ordinary calculus, Eq. fIDll) is not well-defined be- 
cause r{t) is discontinuous at every point. One way to circumvent this problem is to use 
discrete time steps of infinitesimal length 6t rather than continuous time. Then, Eq. flDip 
takes the following form: 

^«^±§^ = .W.(«)| + 9hH0ir(0. (D2) 

However, this is not the only way to interpret Eq. fIDll) . For example, an alternative way to 
go from the continuous to a discrete description is to write Eq. (IDip as. 



r — cu 

6t 



ipit) +ij{t + St) 



+ 9 



ij{t) + ^{t + 5t) 



T{t). (D3) 



2 

In fact, there are an infinite number of ways to interpret Eq. (IDll) . depending on the relative 
weight of ip{t) and ip{t + 6t) inside the arguments of the functions on the right hand side 
of the equation. The two most commonly used interpretations are Ito's and Stratonovich's 
prescriptions. The former corresponds to Eq. (]D2p . and the latter to Eq. (]D3|) . 

In physics, Stratonovich's prescription is commonly used because r{t) is usually an ap- 
proximation to a thermal noise with small but finite correlation time; therefore, the argument 
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of g{-) should be an average value of i/) over the time that the correlations persist. In popula- 
tion genetics, on the other hand, Ito's prescription is appropriate because a random change 
of the allele frequencies depends only on the genetic composition of the population prior to 
the change. 

Without the stochastic term, Eqs. (]D2p and (1D3P would yield the same results provided 6t 
is sufficiently small, but the stochastic terms remain different even in the limit 6t 0. 
An easy way to see this difference is to average Eqs. flD2|) and flD3l) with respect to the 
nondifferentiable noise function r(t). Ito's prescription gives — = {u![ip{t)])St 

because {g[4'(t)]T(t)) = {g[ip(t)]){T(t)) = due to the independence of ipit) and T(t). A 
similar simplification, however, cannot be applied to Stratonovich's prescription because, 
generically, the stochastic term depends on ipit + 6t), which is not independent of r(t). 

Because of the aforementioned ambiguity of the stochastic differential equations with 
multiplicative noise in the continuous time limit, care must be taken while differentiating 
stochastic variables. While the rules of ordinary calculus apply to Stratonovich's prescrip- 
tion, special rules of the Ito calculus are required for Ito's prescription when tracking the 
evolution of, say, a polynomial function u[ip(f}\ of the stochastic variable obeying Eq. ( IDll) . 



In this paper, we use Ito's formula, namely 
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22 



where m(V') is a twice continuously differentiable function, and the primes now indicate 
differentiation with respect to ip. 



APPENDIX E: A MODEL WITH SEVERAL NEUTRAL ALLELES 

A model with q neutral alleles is an intermediate case between the two-alleles model 
that we focus on in this paper and the infinite alleles model discussed in Appendix O The 
g-alleles model is also analogous to nonequilibrium g-state Potts models. In this appendix, 
we briefiy outline how the g-alleles model can be formulated and solved in the language of 
one and two-point correlation functions, compare our analytical predictions to simulations, 
and extend Eq. fH3l) to the undulating front model. 

To specify the g-alleles model, we let fi{t,x) be the frequency of allele i at time t and 
position x; these quantities satisfy 'Y^l^i fiit^x) = 1. The spatial diffusion and coalescence 
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probability of two lineages are still characterized by Dg and Dg respectively. Intra-allelic 
mutations are described by the mutation matrix Hij, which is the probability of allele i 
mutating into allele j Hi 7^ j. When i = j, we let fiu = — X]j=i j^it^ij describe the 
outflow of alleles from allelic state i due to mutations. 

The dynamics of the g-alleles model can be analyzed by considering one-point corre- 
lation functions Fi{t,x) = {fi{t,x)) and two-point correlation functions Fij{t, Xi,X2) = 
{fi{t,xi)fj{t,X2))- Fi{t,x) is the probability to flnd allele i at position x at time t, 
and Fij{t, Xi, X2) is the probability to simultaneously find at time t allele i at position xi and 
allele j at position X2- The evolution equations for these correlation functions are obtained 
by tracing one and two lineages backward in time; the results are 

'-^-n.'^,±,M (El) 

dF,j{t, xi,X2) ^_d_ _^ Fij{t, xi,X2) + Dg5{xi - X2)[SijFi{t,xi) - Fij{t, xi,X2)] + 
q Q 

[fj'i'iFi'j (t, xi, + fj'j'jFij' (t, xi, X2)], 

i'=i j'=i 

(E2) 

where 6ij is Kronecker's delta, which is zero if z 7^ j and one otherwise. Thus, for a generic 
mutation matrix fiij one has to solve a system of coupled linear partial differential equations. 

For simplicity and the ease of comparison with the other results in this paper, let us 
assume spatial homogeneity and identical mutation rates between any two alleles, fii^j = 
11/ q. Under these assumptions, Eq. (]E2[) can be simplified by introducing averaged spatial 
heterozygosity 

i=l j=l i=l 

which is the probability to sample two different alleles at time t distance x apart. The 
equation of motion for H{t,x) can be derived both from Eq. flE2p and, more simply, by 
tracing two lineages backward in time: 

dl^ = ^^^9^2^ + 2^ -Hj- DgH6{x). (E4) 
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Note that Eq. (lE4p agrees with Eq. (1C1[) in the hmit q oo and with Eq. (H6l) for /ii2 = 
Ai2i = yu/2. Since Eq. flE4p has the same functional form as Eq. fICip . the methods of 
Appendix [Q can be used to solve for H(t,x). 

In the absence of mutations, Eq. flE4p is identical to Eq. fl5^ . as we briefly mentioned in 
Sec. IIVI However, g-alleles models with different q may have slightly different dynamics due 
to g-dependent initial conditions: for example, an initially well-mixed population is repre- 
sented by H{0, x) = Hq = 1 — 1/q. Thus, the results of Sec. [IV] apply to the g-alleles model, 
provided appropriate initial conditions are used. In particular, we expect the standard de- 
viation of fi(t), the total frequency of allele i, in a finite population to grow as t^^^. This is 
indeed confirmed by our simulations shown in Fig. (TH Spatial correlations in the nonequi- 



librium g-state Potts model have recently been analyzed by Masser and ben-Avraham |34| . 
who also found that two-point correlation functions obey the same g-independent equation 
of motion. 

Finally, one can obtain the behavior of the standard deviation of the total frequency 
of allele one, A(t), in the undulating front model by the following scaling argument. We 
consider A(t) at large times after monoallelic domains have formed. Let Nd{t) be the number 
of domains consisting of allele one and dk{t), k = 1,2, Nfi{t) be lengths of these domains. 
Then, A(t) is given by 



A(t) 



Nd{t) Na(t) 

5^4(t)-($^4(t)) 

k=l k=l 



(E5) 



We simplify Eq. flESp by making an approximation that N^it) and dkif) for /c = 1, 2, N^if) 
are independent random variables, which gives 

m ^ {{N,m{di{t)) - {d.m + {d.imiNUt)) - {N.m} , (eg) 

where we used the fact that di{t) are identically distributed. 

By using first passage time analysis discussed in Ref. [35|, one can show that {Nj(t)) — 
{N,{t))' oc {N,{t))[{dl{t)) - {d,{t)f]/{d,{t)f oc miit)) - {d,{m/{d,{t)y. Thus, 

AW«;^^[(rf?W)-(^iW)']- (E7) 
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FIG. 14: (Color online) Genetic drift in a finite population with three alleles, (a) The genetic 
composition of the population [fi(t), f2(i), fsCO] projected on the plane Yl^=iUi't) = 1 in a single 
run of the neutral 3-alleles model with a flat front. Here L = 1000, and there are no mutations, 
(b) The standard deviation of the frequency of allele one A(t), shown with blue dots, is obtained 
from 110 realizations of the simulations described in a. The red solid line shows the slope expected 
from Eq. (f43]) . 



Upon recalling, that, in the undulating front model, {di{t)) oc t^, and {di{t)y oc t^^, we 
conclude that 



A(t) oc — oc . 



(E8) 



where, in the last proportionality, we used ( = 2/3 from Ref. [ll| 
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